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1 INTRODUCTION 



^* The Lya forest (Lynds 1971, Sargent et al. 1980, see Rauch 
. {_( | 1998 for a review) arises naturally in cosmological structure 
■ formation scenarios where gravitational instability acts on 
' small initial density perturbations. In hydrodynamic sim- 
5^ i ulations of such models (Cen et al. 1994, Zhang et al. 
1995, Hernquist et al. 1996, Wadsley & Bond 1996, The- 
uns et al. 1998, see also the analytical modelling of e.g., Bi 
1993, Reisenegger and Miralda-Escude 1995), most of the 
absorption seen in high redshift QSO spectra is generated by 
residual neutral hydrogen in a continuous fluctuating pho- 
toionized intergalactic medium. In such a picture, absorbing 
structures have a large physical extent. Observational sup- 
port for this has come from comparison of the Lya forest 
in adjacent QSO lines of sight (Bechtold et al. 1994, Din- 
shaw et al. 1994, 1995, Crotts & Fang 1998). For matter 
in this phase, it is predicted and found in simulations that 
the underlying mass density field at a particular point can 
be related to the optical depth for Lya absorption (see e.g., 
Croft et al. 1997) and hence a directly observable quantity, 
the transmitted flux in the QSO spectrum. 

Much work has been devoted to studying the statistical 



ABSTRACT 

We present predictions for the one-point probability distribution and cumulants 
of the transmitted QSO flux in the high redshift Lyman-ev Forest. We make use of 
the correlation between the Lyman-a optical depth and the underlying matter den- 
sity predicted by gravitational instability theory and seen in numerical hydrodynamic 
simulations. We have modelled the growth of matter fluctuations using the non-linear 
shear-free dynamics, an approximation which reproduces well the results of perturba- 
tion theory for the cumulants in the linear and weakly non-linear clustering regime. 
As high matter overdensities tend to saturate in spectra, the statistics of the flux dis- 
tribution are dominated by weakly non-linear overdensities. As a result, our analytic 
approach can produce accurate predictions, when tested against N-body simulation 
results, even when the underlying matter field has rms fluctuations larger than unity. 
Our treatment can be applied to either Gaussian or non-Gaussian initial conditions. 
Here we concentrate on the former case, but also include a study of a specific non- 
Gaussian model. We discuss how the methods and predictions we present can be used 
as a tool to study the generic clustering properties of the Lya forest at high-redshift. 
With such an approach, rather than concentrating on simulating specific cosmological 
models, we may be in the position to directly test our assumptions for the Gaussian 
nature of the initial conditions, and the gravitational instability origin of structure 
itself. In a separate paper we present results for two-point statistics. 

properties of the mass density field and the generic predic- 
tions of the gravitational instability picture. With the Lya 
forest as a probe of the density, we avoid many of the un- 
certainties associated with the use of the galaxy distribution 
to test theories. In principle, it should be possible, by com- 
bining our theoretical knowledge of gravitational clustering 
with observations of Lya absorption, to test the Gaussianity 
of the initial density field, the picture of Lya formation, and 
the gravitational instability scenario itself. In this paper, 
we will concentrate on one-point statistics, namely the one 
point probability distribution function (PDF) of the trans- 
mitted flux, and its moments. We will use as a tool the 
spherical collapse or shear-free model for the evolution of 
density perturbations which Fosalba & Gaztanaga (1998a 
hereafter FG98, 1998b) have shown to be a good approxi- 
mation to the growth of clustering in the weakly non-linear 
regime and which we find works well in the density regime 
appropriate to the study of the Lya forest. Two point statis- 
tics, which probe the scale dependence of clustering will be 
examined in an accompanying paper ( Gaztanaga & Croft 
1999, Paper II). 



From observations of the Lya forest, we can measure 
the PDF of the flux and its moments. The high resolution 
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spectra of the forest taken by the Keck telescope (Hu et al. 
1995, Lu et al. 1996, Rauch et al. 1997, Kirkman & Tytler 
1997, Kim et al. 1997) allow us to resolve structure in the 
flux distribution, and make high precision, shot noise-free 
measurements of these flux statistics. Here we will use the 
statistical properties of the matter distribution, p(x), to pre- 
dict these observable quantities. 

There are a number of studies which predict the evolu- 
tion of the clustering of density fluctuations, and in partic- 
ular of the PDF. The Zel'dovich Approximation (ZA) was 
used by Kofman et al. (1994). Althought the ZA reproduces 
important aspects of non-linear dynamics, it only results in 
a poor approximation to the PDF and its moments. This 
can be quantified by noticing, for example, that the hierar- 
chical skewness S3 = £ 3 /£ 2 in the ZA is S3 — 4 (at leading 
order in £) instead of the Perturbation Theory (PT) result 
S3 = 34/7 (see e.g., Peebles 1980). One way to improve on 
this is to use the PT cumulants to derive the PDF from the 
Edgeworth expansion (Juszkiewicz et al. 1995, Bernardeau 
& Kofman 1995) . In this case the PDF is predicted to an ac- 
curacy given by order of the cumulants involved. Protogeros 
& Scherrer (1997) introduced the use of a local Lagrangian 
mapping (that relates the initial and evolved fluctuation) as 
a generic way to predict the PDF. In this case, the PDF 
is obtained simply by applying a change of variables (the 
mapping) to the PDF of the initial conditions. The best 
of these two approaches is obtained when the Lagrangian 
mapping is taken to be that of spherical collapse (FG98), 
which recovers the PT cumulants to arbitrary order in the 
weakly non-linear regime. There is yet another possibility, 
which involves performing a perturbative expansion and di- 
rectly relating the moments of the flux to the moments of 
the mass (along the lines proposed in a different context 
by Fry & Gaztanaga 1993). This approach does not use the 
density PDF, and could incorporate more exact calculations 
for the (non-linear) density moments. 

Our plan for this paper is as follows. In Section 2 we 
outline the physical basis for the relation we adopt between 
Lya optical depth and the mass distribution. In Section 3 we 
describe our model for following the evolution of the PDF of 
the density and flux, using non-linear mapping relations pre- 
sented in appendix Al. The cumulants of the flux distribu- 
tion predicted by fully non-linear dynamics are described in 
Section 4, together with the predictions of perturbation the- 
ory. The modelling of the effects of redshift distortions and 
thermal broadening is also described. In Section 5, we com- 
pare our analytical results to those measured from simulated 
spectra, generated using N-body simulations. In Section 6, 
we discuss the effects of non-Gaussian initial conditions, the 
redshift evolution of the one-point flux statistics, and the 
bias between flux and mass fluctuations. We also compare 
to other work on the statistics of the Lya forest flux. Our 
summary and conclusions form Section 7. 



2 LYMAN-ALPHA ABSORPTION AND ITS 
RELATION TO THE MASS DISTRIBUTION 

As mentioned in Section 1, the model we use to relate Lya 
absorption to the distribution of mass is motivated by the re- 
sults of numerical simulations which solve the full equations 
of hydrodynamics and gravity, some including star forma- 



tion in high density regions. It was found in these simula- 
tions (e.g., Hernquist et al. 1996) that most of the volume of 
the Universe at high redshift [z > 2, see Dave et al. 1999 for 
the situation at later times) is filled with a warm (10 4 K), 
continuous, gaseous ionized medium. Fluctuations in this in- 
tergalactic medium (IGM) tend to have overdensities within 
a factor of 10 of the cosmic mean and resemble morpholog- 
ically the filaments, walls and voids seen on larger scales 
in the galaxy distribution at lower redshifts. The dominant 
physical processes responsible for the state of this IGM and 
the Lya absorption produced by it were anticipated by semi- 
analytic modelling of the Lya forest (e.g., McGill 1990, Bi, 
Borner & Chu 1992). For completeness, we will summarize 
these processes below. 



2.1 The Fluctuating Gunn-Peterson 
Approximation 

The physical state of most of the volume of the baryonic 
IGM is governed by the photoionization heating of the UV 
radiation background, and the adiabatic cooling caused by 
the expansion of the Universe. The competition between 
these two processes drives gas elements towards a tight re- 
lation between temperature and density, so that 

T = T oP ?(x), (1) 

where p&(x) is the density of baryonic gas in units of the 
cosmic mean. This relation holds well in simulations for 
Pb J$ 10 (see e.g., Katz, Weinberg & Hernquist 1996). Hui & 
Gnedin (1997) have explored the relation semi-analytically 
by considering the evolution of individual gas elements in the 
Zel'dovich Approximation. They find that the value of the 
parameters in equation ([j]) depend on the history of reion- 
ization and the spectral shape of the radiation background, 
and should lie in the narrow range 4000 K < T < 15, 000 K 
and 0.3 < a < 0.6. 

The optical depth for Lya absorption, r is proportional 
to the density of neutral hydrogen (Gunn & Peterson 1965). 
In our case, this is equal to the gas density pb multiplied by 
a recombination rate which is proportional to pbT~ ' 7 . By 
using equation ([!]), we find that the optical depth is a power 
law function of the local density: 

r(x) = Ap b (xf, (2) 

where x is a distance along one axis, taken to the line-of- 
sight towards the QSO (we are working in real-space for 
the moment). Because this result is simply a generalisation 
of Gunn-Peterson absorption for a non-uniform medium, it 
has been dubbed the Fluctuating Gunn-Peterson Approxi- 
mation (FGPA, see Rauch et al. 1997, Croft et al. 1998a, 
Weinberg et al. 1998a). The FGPA amplitude, A, is depen- 
dent on cosmology and the state of the gas so that (e.g., 
Croft et al. 1999), 




Here T is the 

photoionization rate, h — Ho/WO km s _1 Mpc" 1 , and fit, 
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is the ratio of the baryon density to the critical density. The 
FGPA slope, /3 = 2- 0.7a ~ 1.6. 

The FGPA has been tested in simulations (Croft et al. 
1997, Weinberg 1999), and the predicted tight correlation 
found to hold well. The analysis in this paper will involve 
using equation (^) to relate the optical depth to the under- 
lying real-space mass density. We will make predictions for 
the observable quantity, transmitted flux in a QSO spec- 
trum, which we label d>: 



4>{x) 



-Ap b (xf 



(4) 



Equation (^) can be thought of as "local biasing relation" 
between the flux and mass distributions, ft can be seen that 
in this relation, the only spatially varying quantity is Pb(x) 
(ignoring global redshift evolution and assuming a smooth 
ionizing background). Given that the physical processes in- 
cluded in the derivation of the FGPA relation are the dom- 
inant ones, then the clustering properties of the Lya forest 
should be determined mainly by the statistics of ph. The 
emphasis in this paper is therefore on applying our knowl- 
edge of the behaviour of density perturbations to the Lya 
forest. We will use analytical results for the non- linear evo- 
lution of density perturbations in an effort to understand 
the origin of the values of Lya forest observables. The ul- 
timate aim is that with this understanding, measurements 
made from observational data can be used to directly test 
both the gravitional instability hypothesis, and the picture 
of the Lya forest outlined above, as well as throwing light 
on the nature of the primordial density fluctuations. 

An alternative to the approach we adopt here would be 
to use the local relation EqQ to directly reconstruct the 
density field, rather than to predict its cumulants or the 
cumulants of the flux. This reconstructed field could then 
be used to estimate the statistical properties of the density 
(e.g., cumulants) in a straightforward way. This is not how- 
ever simple to do in practice because of the saturation of 
flux in high density regions. Although large changes in high 
density regions have little effect on the statistics of the flux 
(i.e. the cumulants), they will totally change the statistics 
of the density. Any reconstruction technique will therefore 
have to deal with this missing information somehow. One 
approach for dealing with this problem has been presented 
by Nusser & Haehnelt (1998). In the present paper, we make 
use of the important fact that the power-law and exponen- 
tial weighting of the density in the FGPA relation results 
in a flux distribution whose statistical properties are dom- 
inated by the small fluctuations, i.e. the linear or weakly 
non- linear regime. 



2.2 Additional complications 

There are a number of assumptions concerning the relation- 
sip between flux and mass which we should discuss before 
proceeding. First, the above equations apply to the density 
of gas, pb rather than the total density of matter, p, which 
will be dominated by a dark matter component in the models 
we are considering. At the relatively low densities of inter- 
est here, the distribution of gas in simulations does however 
trace the dark matter well. Pressure forces on gas elements 
tend to be small compared to the gravitational forces, and 
non-hydrodynamical N-body simulations can be used to pro- 



duce very similar spectra to the simulations which include 
these pressure effects (Weinberg et al. 1999). Simulations 
do have finite resolution limitations, though, and cluster- 
ing in a dissipationless dark matter distribution with power 
extending to small scales cannot be followed with infinite 
resolution. The N-body only calculations so far used (e.g., 
Croft et al. 1998a) have a resolution comparable to the small 
scale smoothing produced by pressure effects. Hydrodynam- 
ical simulations at high resolution (e.g., Bryan et al. 1998) 
can be used to study this smoothing. In the case of analytic 
work, one can first consider the linear regime. In this case, 
the power spectrum of fluctuations in the gas density, P g (k) 
is a smoothed version of the dark matter power spectrum, 
PuM.(k), so that 



P 3 (k) = 



[i + (k/kj 



(•>) 



where kj is the Jean's wavenumber (see e.g., Peebles 1993). 
In tests of this result, Gnedin & Hui (1998) have shown that 
after reionization, the effective smoothing length is gener- 
ally smaller, and modelling with a different (Gaussian) filter 
tends to give better results when compared with simula- 
tions. The situation in the non-linear regime will be more 
complicated. The Jeans length scales as (pi,) -0 ' 5 , but due to 
the temperature density relation of equation denser re- 
gions also tend to have higher temperatures, more thermal 
pressure, and more smoothing, so that the overall density 
dependence of the Jean's length should be weak. Gnedin 
& Hui (1998) show that filtering the initial conditions of a 
dissipationless simulation with a single scale gives reason- 
able results compared to the full hydrodynamic case (al- 
though worse than their "Hydro-PM" technique, which in- 
volves adding a pressure term to the dissipationless simula- 
tion calculations). In our case, the analytic approximation 
for gravitational collapse which we use allows for filtering 
the evolved density with a top-hat filter in real-space (see 
the next sections and Appendix A1.2). It may be possible 
to vary the smoothing length as a function of density, but 
for reasons of simplicity we use a constant smoothing ra- 
dius for now. Another possibility for the future might be 
self-consistent modelling of the hydrodynamic effects when 
following the evolution of density perturbations. This has 
been done numerically in ID simulations of spherical col- 
lapse by Haiman et al. (1996). 

Second, the FGPA itself will break down in regions of 
high density, because of shock heating of gas, collisional ion- 
ization, star formation, and other processes. We can quantify 
this by appealing to the results of hydrodynamic simula- 
tions. As stated above, the relation has been directly tested 
by Croft et al. (1997), who find that it works well at high 
redshifts, z > 2, on a point by point basis, for pt < 10. When 
we consider statistics that we might want to measure from 
the flux distribution, the situation is even better. For exam- 
ple, we can see using the numbers given for A above at z = 3 
and equation (Q), that optical depth will saturate (cf> < 0.05) 
for p > 3. The physical processes occuring in regions with 
p > 3 are therefore not likely to directly affect what we can 
measure. Of course, there will be indirect effects, for exam- 
ple, supernova winds may inhomogeneeously heat the IGM 
out into lower density regions. Also, the reionization of Hell, 
which is expected to occur around z ~ 3, may cause inho- 
mogeneous heating if it is patchy enough (Miralda-Escude 
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and Rees 1994). Although we expect the volume occupied 
by regions which do follow the FGPA to be overwhelming in 
the high-z Universe, the statistical properties of the absorp- 
tion predicted by analytical gravitational instability theory 
should be useful in testing the validity of this assumption. 
They should also help us decide if their is any apprecia- 
ble contribution to clustering from spatial variations in the 
photoionization rate F, due to the inhomogeneity in the UV 
background. Any such variations are expected to be small 
in amplitude and to occur only on scales larger than we can 
probe directly at present (see e.g., Zuo 1992, Fardal & Shull 
1993, Croft et al. 1999). 

Third, we have so far only dealt with the density field 
in real-space, whereas measurements from QSO spectra are 
made in redshift-space. Both peculiar velocities and thermal 
broadening of absorption features should affect the statistics 
of <j) to some degree. We will include both these effects in 
our predictions. 



3 THE PROBABILITY DISTRIBUTION 

We would like to make predictions for the one-point PDF of 
the flux cj> and its moments. The one-point PDF of a given 
field <f> is defined so that the probability of finding, at a 
random position x, a value <f>(x) in the infinitesimal interval 
(j) to (j> + d(f>, is P(<f))d<j) 

To make these predictions we will first derive the cor- 
responding probabilities P(p), for the local matter overden- 
sity p = 1 + 5, where p is in units of the mean density 
p(x) = n(x)/{n). We will make indiscriminate use of either 
8{x) or p(x) as variables when describing density fluctua- 
tions. The second step is the assumption of a local relation 
with the form (f> = f(p), (motivated by equation jij). The 
PDF of the flux will then simply be obtained by performing 
a change of variables from the PDF of the density. 

We start by assuming a (Gaussian) form for the PDF of 
the initial conditions, and then follow its evolution . As we 
will see, in our approach it is not necessary to assume Gaus- 
sian initial conditions, and this procedure can be extended 
to some other non-Gaussian models. We will do this with a 
model starting from \ 2 initial conditions in Section 6.1. 

One important point to note about our predictions is 
that we are not creating artificial spectra but instead us- 
ing an analytical model to evolve the density PDF and then 
predict the PDF of the flux directly. Our predictions for the 
density distribution will depend on only two parameters, 
equivalent to the slope and amplitude of the linear correla- 
tion function on the smoothing scale, 7, and erf, (see Sec- 
tions 3, 4 and Appendix Al). This will allow us to cover a 
wide range of possiblities, and make predictions that are as 
generic as possible. 

In this section we will also test the effects of varying 
the two parameters in the FGPA relation, A and /3, which 
(equation |^]) contain information about the cosmic baryon 
density, ionizing background and reionization history of the 
Universe. 



3.1 The PDF of the initial conditions 

In the limit of early times, we assume an nearly homoge- 
neous distribution with very small fluctuations (or seeds), 



with given statistical properties. We will concentrate on the 
case where the statistics of the initial density are well de- 
scribed by Gaussian initial conditions, which correspond to 
a general class of models for the initial conditions. The one- 
point Gaussian probability distribution of an initial field 8 
is given by: 



Pic iS) = 



exp -- 



(6) 



As the overdensity must be positive, p > 0, we have 
that 8 > — 1, and a Gaussian PDF only makes physical 
sense when the initial variance is small: ctq — * 0. 



3.2 The evolved mass PDF 

Because of gravitational growth, the evolution of 8 will 
change the PDF from its initial form. For small fluctuations 
linear theory provides a simple way of predicting the time 
evolution of 5(t,x) — D(t)So(x), where D(t) is the growth 
factor (equal to the scale factor D — a for £1 = 1), and So(x) 
is the initial field. We will denote this linear prediction by 
5l- For Gaussian initial conditions the linear PDF is also 
Gaussian with a variance a\, given by scaling the initial 
variance a 2 by D 2 , so that erf = D 2 a 2 . 

As mentioned in the introduction, there are a number 
of studies which predict the evolution of the PDF beyond 
linear theory. Here we will consider a generic class of lo- 
cal mappings along the lines introduced by Protogeros & 
Scherrer (1997). The idea is for us to relate the non-linear 
fluctuation S(q) (in Lagrangian space) with the correspond- 
ing linear fluctuation Sl{<}) = DSo{q) using a universal (lo- 
cal) function. To simplify notation we choose to express this 
mapping as a relation between the non-linear overdensity 
p — 8 + 1 and the linear fluctuation Sl , so that 



p(q) = G[S L (q)]. 



(7) 



One such mapping is the spherical collapse model (SC) 
or shear-free approximation. For Gaussian IC, the SC ap- 
proximation happens to give the exact statistical properties 
of the density (cumulants of arbitrary order) at leading order 
in perturbation theory (as found by Bernardeau 1992 in the 
context of the cumulant generating function), and provides a 
very good approximation to higher orders (see FG98) . Phys- 
ically, this mapping corresponds to taking the limit where 
shear is neglected. In this case the equations for the growth 
of 8, in Lagrangian space, are identical to those given by 
spherical collapse. So, in the perturbative regime, the SC is 
the best mapping possible, given the local assumption made 
in equation (pj). The local transformation naturally occurs 
in Lagrangian space q (comoving with the initial fluid ele- 
ment). The important point to notice here is that although 
the local mapping is not the exact solution to the evolution 
of 8 (which is in general non-local), it does give the correct 
clustering properties in the weakly non-linear regime. In the 
Appendix we give some specific examples for the transfor- 
mation Q. 

The one-point PDF induced by the above transforma- 
tion, in terms of the initial one-point PDF Pic, is 



Pl(p) = Pic {8 1 



d8 L 



dp 



(8) 
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where Sl = <7~ 1 [p]- As mentioned before, the above ex- 
pression corresponds to the probability distribution of the 
evolved field in Lagrangian space, q. To relate Lagrangian 
and Eulerian probabilities we use the law of mass conserva- 
tion: d8(q) = p dS(x), where p(x) = 1 + S(x) is the overden- 
sity in Eulerian coordinates. We therefore have 



P(p) = 



1 Pic(Sl) 
N p 



dS L 



dp 



(9) 



where A r is a normalization constant. We will show some of 
these predictions in Section [| (e.g., Fig. Q). 



3.3 The PDF of the flux 

We next assume a local transformation which relates the un- 
derlying smoothed overdensity to some observable quantity 



= f(p) 



(10) 



This quantity can further be related to the linear density 
field, so that 



4> = f[G{6z.)]. 



(11) 



The PDF of (j> will then related to that of the density 
by a simple change of variable: 



P(<t>) = P{p) 



dp 



1 Pic(Sl) 
N p 



d8 L 



(12) 



where Sl = Q~ \f~ [<f>]] and p = / _1 [<£]. Thus, given the 
transformations / and Q, the above equations provide us 
with analytical (or maybe numerical) expressions for the 
PDF of 4>. 

In the present work, we concentrate on the cumulants 
(see Section 5.3 for a definition) of the flux, rather than 
the PDF itself. The reason for this is that, given the local 
assumption, the cumulants are more accurately determined 
(seeSectionUfor more details). Nevertheless, we will see in 
Section ^ that the above prediction gives a good qualitative 
description of the PDF (see e.g., Fig. H). 



3.4 Redshift-space distortions 

The smoothed flux and its corresponding optical depth 
t = — In (j> has been assumed to be a local function of the 
smoothed non-linear density p. The optical depth r at a 
given real-space position along the line of sight r will lie 
at a redshift-space position s in a QSO spectrum: 



s = r + V r /H, 



(13) 



where v r is the component of the smoothed peculiar velocity 
along the line of sight at r. Note that the redshift distortion 
is of the smoothed field, where the smoothing, due to finite 
gas pressure, occurs in real-space. The redshift mapping will 
conserve optical depth, r s ds = rdr, so that we have: 



1 + 



1 dv r 

H~d7 



(14) 



In general, the relation between dv r /dr and p will be com- 
plex. However, in the SC model, spherical symmetry leads 
to a great simplification: 



dv 



dr 3 V ' 



H 
— i 

3 



(15) 



as, by symmetry, derivatives are the same in all directions 
(this idea has also been used by FG98 and by Scherrer & 
Gaztahaga 1998). We can now again use the local mapping 
to relate velocity divergence to the linear field: 8 = Q v \Sjj\, 
as in equation (k.7). Thus we have that the redshift optical 
depth is given by a different mapping: 



1 



Ts(p) = r(p) 
The redshift-space flux is simply: 
4> s = exp[-T s (p)], 



(16) 



(17) 



and its PDF can be computed with a simple change of vari- 
ables: 



.) = 



1 Pic(Sl) 
N p 



dSr 



(18) 



4 THE CUMULANTS 
4.1 Definitions 

Consider a generic field, <f>, which could be either the mea- 
sured flux in a ID spectrum, cj> = 4>{p), or the mass density in 
3D space (f> — p. The Jth-order (one-point) moments of this 
field are defined by (note the subscript "c" for connected): 



mj 



J ) = / p(4>) 4> J d<t>. 



(19) 



Given the above relations we can choose to calculate the 
moments by integrating over Sl'- 

Pic[8l] 



mj 



d5r 



p{5l 



■<f> J (p(s L )), 



or over the non- linear overdensity p: 



mj 



dp 



Pic[5l(p)] 



P 



d8 L 



dp 



J (P)- 



(20) 



(21) 



Here <p can refer either to real-space fields or fields which 
have been distorted into redshift-space (e.g., equation ]l7|]). 
The Jth order reduced one-point moments, or cumulants kj, 
of the field <f) are defined by: 



kj 



aiog M(t) 



dt 



(22) 



where M(t) = (exp(^t)) is the generating functional of the 
(unreduced) moments: 



mj 



j dM(t) 
' dt 



The first reduced moments are: 



fci 
k 2 
k 3 



nil 

m 3 
m4 



(23) 



(24) 



3mim2 + 2m! 

6mi + 12m\m2 — 3m2 



4mim3 



and so on. Note that even when we normalise the flux so 
that the mean is zero (mi = 0), the cumulants are different 
from the central moments in that the lower order moments 
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have been subtracted from them, so that &4 = 7714 — 3m| for 
mi = 0. 

It is interesting to define the following one-point hier- 
archical constants: 



J > 2 



(25) 



This quantities turn out to be roughtly constant under grav- 
itational evolution from Gaussian initial conditions (see e.g., 
Gaztanaga & Baugh 1995 and references therein). 



4.2 Fully non-linear predictions 

We now use the one-point flux PDF obtained from equa- 
tion (|l^) to predict the one-point moments of the flux. As 
mentioned in Section 3, to make the predictions we need 
the value of the linear variance a\, and the local relation 
equation (Q), which only depends on the smoothing slope 
7 (see equation A17 for the definition of 7). For non-linear 
mapping relations, we try each of the two cases introduced 
in the Appendix. In the following figures we will only plot 
results for the Spherical Collapse (SC) mapping because 
they coincide perfectly with the results for the Gener alized 
Zel'dovich Approximation (GZA) model in equation (A10) 
with a = 21/13. 

We are interested primarily in the evolution of the den- 
sity PDF and the weighting which the FGPA relation gives 
to the clustering properties of the density field. In order to 
separate the effects of redshift distortions from those of den- 
sity evolution, we will present results in real-space first. 



4-2.1 The mean flux 

Figure ^ shows the mean flux {4>} (= mi in equation jti|) as 
a function of a\ for several values of A, f3 (the parameters 
in the FGPA relation) and 7. 

For small o~l all predictions tend to {<fi) — exp(-A) 
which corresponds to the flux at the mean overdensity: <j> — + 
exp(-A) as p — > 1. We will see in next subsection that this is 
the leading order PT prediction. For larger ox the mean flux 
becomes larger, as expected, but it is flatter as a function of 
cr|_ when there is less smoothing, i.e., when 7 is less negative. 
Less smoothing of the density will correspond to larger non- 
linearities, at least in PT (see e.g., FG98). This seems to 
indicate that the effect of non-linearities is to reduce the 
mean flux, competing with linear growth, which increases 
the mean flux. 



4-2.2 The variance of the flux 

We define the variance using the normalized flux: 



(<t>) 



{4>) 



(26) 



The overall normalization by (</>} is a convention, in analogy 
to what is done for density fluctuations. Figure |^ shows the 
predicted variance as a function of a\ for several values 
of A, P and 7. For small a\ these results reproduce the linear 
relation: = b\ a\ (see equation Ji^]). Deviation from this 
power-law relation (of index 1) occurs as a\ is increased, 
and occurs sooner for larger values of A and f3. For 7 > 0, 
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Figure 1. Mean flux (<f>) as a function of the linear variance <r£ 
for different values of A and (3. The dotted, short-dashed and 
long-dashed lines show the predictions for 7 = 0,-1 and —2 re- 
spectively. 
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Figure 2. The variance of the flux o - ^ as a function of <r£. The 
meaning of the line types is the same as in Figure |l|. 



when a\ reaches ~ 1, o"2 seems to reach a maximum and 
then decreases again like a power-law for large a\. 

We can again see that the predictions become flatter as 
a function of o~\ when there is less smoothing, i.e., for less 
negative 7. 



4-2.3 The skewness of the flux 

In a similar way, we define the (normalized hierarchical) 
skewness of the flux as: 
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Figure 3. The Hierarchical skewness of the flux, 53(0), as a 
function of <r£. The meaning of the line types is the same as in 
Figure [l]. 



s 3 (4>) 



(WW 



(27) 



Figure |^ shows the predicted skewness, as a function of erf 
for several values of A, /? and 7. Because of the dependence 
of flux on density in the FGPA relation (more density, less 
flux), while the density distribution is positively skewed, the 
skewness of the flux tends to be negative for most values of 
the parameters. 

For small erf the results tend to a constant as expected 
in leading order PT (see equation ^]). We will examine 
the PT relations in more detail in Section 4.3. For the mo- 
ment, we note that there again seems to be less variation in 
Ss(c6) as a function of erf for cases with less smoothing (less 
negative 7). 

4.2.4 The kurtosis of the flux 

The (normalized hierarchical) kurtosis of the flux is defined 
in a similar fashion: 



(WW - i) 4 ) c 



(28) 



Figure ^ shows the predicted kurtosis, as a function of erf 
for several values of f3 and 7. For clarity we only show one 
value of A. 

For small erf, these results tend to the constant value 
predicted by leading order PT (see equation Q). Being a 
fourth moment, 5*40 is extremely sensitive to deviations from 
Gaussianity which occur as the density field evolves and erf 
increases. This sensitivity seems to be larger for lower values 
of 0, presumably because when /3 is high, high density parts 
of the PDF which might contribute heavily to the kurtosis 
of the density field have little weight in the statistics in <j>. 



Figure 4. Same as Figure ^ for the Hierarchical kurtosis. For 
clarity only A = 0.6 is shown. 



4.3 Perturbative predictions 

An alternative to using the PDF is to calculate the cu- 
mulants directly from the perturbative expansion along the 
lines suggested (in the context of galaxy biasing) by Fry & 
Gaztahaga (1993). That is we take: 



Cfe r fc 



61 

CO 
Cl 
C2 
C3 



-ApP _ 

-Afie 
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~A~P 



2 - 3/3 + f3 2 + 3AP - 3Af3 2 + A 2 f3 2 



(29) 



(30) 



and so on. From this expansion one can simply estimate the 
moments by taking mean values to the powers of <j>. The 
leading order terms in erf are: 

— A 



{4>) 
2 

v.* 



(31) 



,2 2 

53 + 3c 2 

61 

5 4 + 12ff 3 C2 + 12cj + 4c 3 



(32) 



where £3 and 54 are the leading order (hierarchical) skew- 
ness and kurtosis for the density field. For Gaussian initial 
conditions they are: S3 = 34/7 + 7 and S 4 = 60712/1323 + 
62/37 + 7/37 2 (both in the SC model and PT theory). For 
non-Gaussian initial conditions, one would have to add the 
initial contribution, e.g., S3 = S3 + 34/7 + 7. 

These predictions are shown as a straight continuous 
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Figure 5. Perturbative predictions for the one-point flux mo- 
ments compared to the fully non-linear predictions. Mean flux (<j>) 
(top), variance (X^ (middle panel) and skewness S$(<f) (bottom) 
are shown as a function of the linear variance cr^ for A = 0.6 
and (3 = 1.0. The dotted and long-dashed lines show the non- 
linear predictions (Section 4.2) for 7 = and 7 = —2 smoothing 
respectively. The straight continuous lines show the correspond- 
ing leading order perturbative prediction (Section 4.3), valid for 
ox — > 0. The solid curve in the top panel is the perturbative 
prediction for (<f>) including the effect of a higher order (loop) 
correction (see equation [B3j). 

lines in Figure ||, where they are compared to the full (non- 
perturbative) calculation in the SC model for two values of 
the smoothing slope, 7 = (long-dashed line) and 7 = —2 
(dotted line). We can see that the expressions only converge 
on the correct result asymptotically as ox — ■> 0. The relative 
performance depends on 7, with steeper slopes giving better 
results. It is easy to calculate higher order (loop) corrections 
(see e.g., FG98). For example: 



1 + 4^(1-/3 + 49) al + 0(a 4 L ) 



(33) 



This prediction for the mean flux is shown as a curved con- 
tinuous line in the top panel of Figure ||, and seem to work 
up to scales where <tl — 1. We find however that the agree- 
ment becomes worse for larger values of A and /3. A similar- 
tendency is found for other moments. 

We have seen that even when high-order corrections 
are included, this perturbative approach only works well for 
small values of ox , A and 0. It is likely to have only a limited 
applicablity when we consider the situation in the observed 
Universe, where typical values of ox > 1 are expected, (at 
least for redshifts z < 4). Given that we have the possibil- 
ity of implementing the SC model (or the GZA model) to 
arbitrary order, one could ask why we bother with a pertur- 
bative approach. The first obvious reason is that it gives us 
compact analytical predictions, simple formulae which are 
functions of the input variables (A, (3, 7 and cr£). A second 
reason is that by using this approach, it may be possible to 
introduce the PT solutions. As mentioned before, PT only 



Figure 6. The effect of redshift-space distortions on the one- 
point moments. Mean flux (</>) (top), variance ct^ and skewness 
Ss{(j>) (bottom) are plotted as a function of the linear variance 
<T? for 7 = 1, A = 0.8 and /3 = 1.6. The short-dashed and 
continous lines show the predictions in real and redshift-space 
(peculiar velocities only) respectively. The long-dashed lines cor- 
respond to predictions in redshift-space with an additional ve- 
locity dispersion component on small scales (added in order to 
match simulation results - see text). The dot-dashed lines show 
the predictions in redshift-space with thermal broadening as well 
as peculiar velocities (and the extra small-scale dispersion). In 
modelling the thermal broadening component to the redshift dis- 
tortion we have assumed that the temperature depends on the 
density, as predicted by equation (y) (see Section 4.4). The points 
show results from simulated spectra (set [a]) described in Sec- 
tion 5. The triangles, circles and squares represent spectra in real 
space, redshift-space with no thermal broadening, and redshift- 
space with thermal-broadening, respectively. 



differs from the SC model through the shear contributions, 
and although these are typically small (as will be shown) 
they might still be relevant for higher precision comparisons. 
It is not clear nevertheless than one could obtain higher ac- 
curacy at o_l > 1, given the limitations of a perturbative 
approach (i.e., convergence of the series). 



4.4 Predictions in redshift-space 

We now turn to the more observationally realistic case where 
redshift distortions are included. Fig. || shows how the (fully 
non-linear) predictions change when given in redshift-space. 
We use the formalism of Section 3.4 (e.g., equation ]l7]]), 
results which allow one to estimate the effects of peculiar 
velocity distortions. The redshift distortions caused by ther- 
mal broadening can be treated in a similar way, and we will 
also do this below. For clarity we only show a single value 
of A, (3 and 7, but similar effects are found if we use other 
values. We have also plotted some simulation points on Fig. 
|i| The simulations will be decribed fully in the next section. 
For the moment, it is only necessary to mention here that, 
for purposes of comparison, the simulation density PDFs 
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can be decribed by the two parameters, 7 and a\ (evalu- 
ated from the power spectrum of initial fluctuations used 
to set up the simulation). The transformation from density 
into flux in the simulations has also been carried out us- 
ing the FGPA relation. We plot points both with and with- 
out including redshift distortions from peculiar motions and 
thermal broadening. 

If we concentrate on the mean flux first, and ignore ther- 
mal broadening, we can see that {</)) in redshift-space (conti- 
nous line in the top panel) seems to converge to same value 
as the real-space mean flux (short-dashed line) for small ol- 
For <Ti > 1 the mean flux is larger in redshift-space. 



The variance in the flux, a^, defined by equation ( J26[ ) 
is shown in the middle panel of Fig. [| We can see that 
this quantity is larger in redshift-space than in real-space 
for <7i < 1, a trend that is reversed for o\l > 1. The former 
is presumably due to the same "squashing effect" that is 
seen in studies of the density field, where infall of matter 
into high density regions enhances clustering (Kaiser 1987). 
The latter effect can be attributed to a relative decrease in 
the level of redshift-space clustering caused by high velocity 
dispersion along the line of sight. 

The skewness, defined by equation ( |27| ) is shown in the 
bottom panel of Fig. ^j. For small ox, the redshift-space 
(continous line) S3(<j>) seems to match the real-space value 
(dashed line). For ox ^ 1, S3(4>) is smaller in redshift-space. 

The simulation points on the plot in real space (triangle) 
agree with predictions in the real space. In redshift-space 
(open circles), although the sign of the change caused by 
the peculiar velocity distortions is correct, the predictions do 
not agree in detail. Our interpretation of this is that the SC 
model does not predict enough random non-linear velocity 
dispersion. We have therefore added a velocity dispersion 
by hand, by adding in a velocity divergence term, #dis P to 
equation (117 



T s (p) = r{p) 



l + -[e(p) + e diBp (al)] 



(34) 



In order that the asymptotic behaviour of #dis P be satisfied 
(non- linear dispersion — > as o\ — *• 0), we have adopted 
the functional form #dis P = Ca\, Predictions from the SC 
model including the effect of this term are shown in fig. ^ as 
a long dashed line. We have adjusted the constant C so that 
the predictions go through the simulation point in redshift- 
space. 

Next, we include thermal broadening in our predictions. 
For the moderate optical depths we are concerned with here, 
the relevant Voigt profile can well be approximated by a 
Gaussian velocity dispersion. The width of this Gaussian 
profile, o T ^ a T[> (T/T0) 1/2 , where ctto lS/v^kms" 1 
for T ~ 10 4 K. From equation (|TJ> we have T oc p 6 , so 
that or — drop ' 3 - We can think of thermal broadening as 
resulting in the addition of a thermal velocity component, 
8t, to the divergence field 9 in equation (|l6|). We will model 
this thermal component in a similar way to the extra non- 
linear dispersion term (#dis P ) defined above, which can be 
thought of as a "turbulent" broadening term. We simply 
model the additional thermal dispersion using its rms value, 
so that, 



<P) 



WK 



O"T0 0.3 



(35) 



where A is distance in the QSO spectrum which corresponds 



to the scale of Jean's smoothing. This density-dependent 
term then enters the RHS of equation (^), alongside #dis P - 
The long-dashed line in Fig. shows the effect of thermal 
broadening using this prescription. We have used the value 
<tto appropriate to the simulation (see Section 5.2 for de- 
tails), whose thermally broadened results are plotted as a 
solid square. As can be seen, thermal broadening results in 
more flux being absorbed, and yields a lower value of (A. 
This is as we should expect, given that the distribution of 
optical depth has effectively been smoothed out by the ad- 
dition of a dispersion. 

It is evident from these results that the one-point mo- 
ments depend fairly strongly on the details of redshift dis- 
tortion modelling, which are likely to be poorly constrained 
a priori in our approximate treatment. Fortunately, much 
of the uncertainty can be removed by setting the mean flux, 
((/>), to be equal to some (observed) value. We will show later 
(Section 5.5) that by doing this, the moments can be made 
insensitive to the inclusion of redshift distortions. 



5 COMPARISON TO SIMULATIONS 

As the emphasis of this paper is on the role of gravitational 
evolution of the density field, we now test the analytic tech- 
niques we have employed against numerical simulations of 
gravitational clustering. The N-body only simulations which 
we use do not allow us to perform tests of the validity of the 
model we have assumed for relating the mass distribution 
to an optical depth distribution (the FGPA, Section 2). Any 
difference between the statistical properties of the Lya for- 
est we predict and those observed in nature could therefore 
stem from a misapplication of the FGPA, rather than from 
problems with the underlying density field. Tests performed 
in other contexts show that this is unlikely, as a dissipation- 
less approach to simulating the Lya forest can perform well 
(see e.g. Croft et al. 1998a, Weinberg et al. 1999) in compari- 
son with the full hydrodynamic case. Approximate methods 
should neverthless be tested on a case by case basis, and 
we reserve comparisons with hydrodynamic simulations for 
future work. 



5.1 Simulations 

The simulations we use have all been run with a P 3 M N- 
body code (Efstathiou & Eastwood 1981, Efstathiou et al. 
1985). The softening length of the PP interaction was made 
large (1 mesh cell), for simulation sets (a-c) because high 
spatial resolution is not needed. We have not attempted 
to simulate any particular favoured cosmological models or 
even to make sure that cases are likely to be compatible 
with expectations for the nature of the density field at high 
redshift. We are more interested in spanning a wide range of 
values of a\, and 7, the parameters which determine the na- 
ture of the evolved density field in the SC model. To this end, 
we use outputs from three different sets of simulations. The 
initial conditions for all runs were Gaussian random fields 
with CDM-type power spectra of the form specified by Ef- 
stathiou, Bond & White (1992), with a shape parameter F, 
so that 

P ( k ) « m , 7ZZ I ,uJ„ , 7^WT„ (36) 



[1 + (ak + (Wc) 3 / 2 + (cfc) 2 )"] 2 /" 
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where v = 1.13, a = 6.4/T 7i _1 Mpc, 6 = 3.0/r ^ _1 Mpc, 
c = 1.7/r /i _1 Mpc. There are five realizations with differ- 
ent random phases in each set of simulations, which are de- 
scribed below. 

(a) A set with a box-size 40 /i _1 Mpc and shape parame- 
ter F = 0.5. The model was run to z = 3 with an Einstein-de- 
Sitter cosmology. At that redshift a\ — 2.0. at the smooth- 
ing radius (see below), which was 0.31 /i -1 Mpc (comoving). 
The linear slope at the smoothing scale, 7 = —0.8. These 
simulations were run with 200 3 particles and 256 3 cells. 

(b) A set with box size 22.22 /i _1 Mpc, 128 3 particles, 
r = 0.5, and an Einstein-de-Sitter cosmology. Simulated 
spectra were made from this set assuming that z — 3, 
but several different outputs were used with varying am- 
plitudes of mass fluctuations. The smoothing radius was 
0.31 h~ 1 Mpc (comoving), on which scale the linear slope 
was 7 = —0.8. The value of erf, on this scale ranged from 
0.02 to 7.5 for the different outputs. 

(c) A set the same as (b), except with T — 10. The 
smoothing radius was again 0.31 /i _1 Mpc (comoving), on 
which scale the linear slope was 7 = —1.8. The value of erf 
on this scale ranged from 0.03 to 10 for the different outputs. 

(d) A set with box size 20 ft~ 1 Mpc, 126 3 , particles and 
r = 3.8. This set was originally used as ensemble G in 
Baugh, Gaztanaga & Efstathiou 1995, but the box size was 
taken to be larger, and hence F was lower. The smoothing 
scale is 0.37 /i~ 1 Mpc, where 7 ~ —1.5, and erf — 0.25 (at 
2 = 3). The cosmology assumed has Q = 0.2 and Qa = 0.8 
at z — 0. 

We have also run a single simulation with the same 
parameters as those in set (a), except with a box of size 
11.11 /i _1 Mpc, and 128 3 particles, so that the mass resolu- 
tion is increased by a factor of ~ 12. 

5.2 Simulated spectra 

To make simulated spectra from the N-body outputs, we use 
the following procedure (see also Hui, Gnedin & Zhang 1997, 
Croft et al. 1998a). The particle density and momentum dis- 
tribution is assigned to a (256 3 ) grid using a TSC (triangular 
shaped clouds, Hockney & Eastwood 1981) scheme. The re- 
sulting fields are smoothed in Fourier space with a filter, 
which in our case is a top-hat with radius given in Section 
5.1 above. We also use a very narrow Gaussian filter (with er 
0.1 times the top-hat radius) in order to ensure that density 
is non-zero everywhere. The velocity fields are computed by 
dividing the momentum by the density everywhere. Spec- 
tra are selected as lines-of-sight through the box, parallel to 
one of the axes (we select the axis randomly for each line-of 
sight). These one-dimensional density fields are then con- 
verted to an optical depth distribution using equation (Q). 
We also compute the temperature distribution of the gas us- 
ing equation (|l|), with a = 0.6 and To = 10 4 K. The optical 
depths are then converted from real-space to redshift-space 
by convolution with the line-of-sight velocity field and with 
a Gaussian filter with the appropriate thermal broadening 
width. 

We estimate the one-point statistics of the flux in the 
spectra without any additional smoothing using counts-in- 
cells and the estimators of Section 4.2. We extract 5000 spec- 
tra from each simulation realization and estimate statistical 
errors from the scatter in results between the 5 realizations. 




Figure 7. The one-point PDF of the density field p measured 
from N-body simulations (continuous line, set [a], see Section 5.1). 
The simulation density field was smoothed with top-hat cell on a 
scale with linear variance erf ~ 2 and slope 7 = — 1. The predic- 
tion of two approximations to PT are also shown: the Spherical 
Collapse model (short-dashed line) and the GZA in a = 21/13 
dimensions (long-dashed line). 

We find that the resulting error bars are typically much 
smaller than the symbol size and so we do not plot them in 
the figures. Small systematic errors will arise because of the 
additional smoothing involved in using a mass assignment 
scheme, also because of shot noise from particle discreteness. 

5.3 The PDF of density and flux 

We first compare the PDF of the density in the simulations 
with the analytical predictions, both in real space. The one- 
point density PDF estimated from simulation lines-of-sight 
(using simulation set [a]) is shown in Fig. |asa continuous 
line. The analytical predictions are shown as a short-dashed 
(SC) and long-dashed (GZA) line. In evaluating the predic- 
tions we have used 7 ~ — 1 and erf ~ 2 which correspond to 
the appropriate linear theory values for simulation set (a). 

Note that the spherical collapse model is only an ap- 
proximation to PT, so we do not expect to recover the PDF 
exactly, even close to 8 ~ 0. To carry out the exact recovery 
we would have to include non-local (tidal) effects. The mean 
tidal effects vary in proportion to the linear variance and the 
leading contribution (when variance goes to zero) is only lo- 
cal (but non-linear). Tidal effects seem to distort the PDF, 
turning some 8 = fluctuations into either voids 8 ~ — a 
or overdensities 8 ~ a, so that the real peak in the PDF is 
slightly lower than SC or GZA. This distortion is significant, 
given that the statistical errors on the simulation result are 
small. The overall shape is however similar, as we can see 
from Fig. ^[ The lower order moments of the density are 
also fairly similar: tidal effects tend to cancel out. When the 
variance is small (and the PDF tends to a Gaussian), the 
PDF can be defined one-to-one by its moments (e.g., with 
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Figure 8. PDF of the flux evaluated using the same density 
distribution plotted in Fig. Three different values of the FGPA 
parameters A and f) are used as shown in each panel. The one- 
point PDF in the simulations is shown as continuous lines and 
is compared with the predictions of the two PT approximations, 
the Spherical Collapse model (short-dashed line), and the GZA 
in a = 21/13 dimensions (long-dashed line). All results are in 
real-space. 



Figure 9. PDF of the flux evaluated from the same density 
distribution used in Fig. |?| (simulation set [a]). The one-point PDF 
in the simulations is shown as a solid line. Real-space results are 
plotted in the bottom panel, and results in redshift-space with 
thermal broadening in the top panel. These results are compared 
with the appropriate analytical predictions, using the Spherical 
Collapse model (short-dashed line) and the GZA in a = 21/13 
dimensions (long-dashed line). 



an Edgeworth expansion). In this limit tidal effects vanish 
and both the PDF and the moments are given exactly by 
the SC model. 

Thus, for the density distribution, the statistics of the 
moments are dominated by the contribution of local dynam- 
ics to the PDF and tidal effects are subdominant (they tend 
to cancel out when taking the mean). It is plausible that a 
similar sort of cancellation will happen for the statistics of 
the flux. In the next sections, we will therefore concentrate 
on predictions for the moments of the flux. 

Different density regimes will be prominent if we con- 
sider the PDF of the flux, as it is a transformed version of 
the density PDF. In Fig. |] we apply the FGPA relation (ecu- 
ation 4) in real-space to the density PDF from Fig. [?]. The 
resulting flux distributions are rather flat, with the high den- 
sity tail being confined to a small region of flux space near 
(f> — 0. Varying the parameters A and (5 produces most no- 
table differences in the fraction of spectra which show little 
absorption. We can see that the lower f3 is, the less likely 
it is that any pixels will be seen with values near the un- 
absorbed QSO flux. The specific case with (3 = TO is not 
realistic, though, as (3 is expected to be > 1.6 for all reason- 
able reionization histories (Hui & Gnedin 1997). 

In Fig. ^ we show the effect of redshift distortions on 
the flux PDF. We have used the prescription of Section 4.4 
for making the predictions, including thermal broadening. 
The effect of the distortions is to evacuate the low density 
regions, as we might expect. The effect is very noticeable, 
as the peak of the PDF is raised substantially. Note that to 
to make this plot, we have not added an additional small 
scale velocity dispersion, as we did previously to obtain a 



closer match to the one-point moments in simulations (Fig. 
^|). We can see that the agreement of the predictions and 
simulations for the PDF in Figure is still good, showing 
that our redshift modelling is at least a good qualitative 
approximation to the underlying physical processes. 

5.4 Cumulants of the flux 

In Fig. ^we compare predictions and simulation results (for 
simulation sets [a] and [d]) for the mean flux {(j>), the vari- 
ance a 2 ,, skewness 53 (0) and kurtosis 54(c/>). In all cases we 
have used /3 = 1.6, and we show the moments as a function 
of the value of A. All results in this fig. have been evalu- 
ated in real-space. Squares correspond to simulation set (a), 
in which the one-point linear variance of the density field, 
erf, = 2 (the non-linear value is a 2 ~ 4) and the linear slope 
on the smoothing scale is 7 = — 1. Triangles correspond 
to set (c), for which a\ — 0.25 and 7 = —1.5. We can see 
that the predictions (from the SC model) are in good overall 
agreement with the simulations. It is encouraging that the 
agreement is not noticeably worse for set (a), which has a 
much higher amplitude of mass fluctuations on the relevant 
scale. 

We can also compare the predictions of the SC model 
with simulations for different values of the amplitude of mass 
fluctuations. In Figures [H] and ^ we have done this, by plot- 
ting the flux moments for several different simulation output 
times. We can see that for most of the range, the predictions 
work well, for both sets of simulations with different linear 
slopes (7). At very low amplitudes erf,, the simulations will 
still be dominated by the effect of the initial particle grid, so 
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Figure 10. The first 4 moments of the flux in simulation set 
(a) (squares) and set (d) (triangles) compared to the appropriate 
analytical predictions (lines), as a function of A. The top panel 
shows the mean flux (cf>) (closed figures) and the variance <r 2 (open 
figures). The bottom panel shows the hierarchical skewness 53(0) 
(closed figures) and the kurtosis S-i (o) (open figures). All results 
are in real-space. 



that the differences we see for the lowest amplitude output 
are not surprising. For high values of a\ , the SC predictions 
start to break down, with S3 and S4 being most affected. 
This shows that we must be careful when intepreting results 
which appear to indicate high values of o\. 

Ultimately, the real test of how far we should trust our 
analytic methods should come from comparisons with hy- 
drodynamic simulations, and in particular those that been 
run at resolutions high enough to resolve the Jean's scale. In 
such simulations (e.g., Bryan et al 1998, Theuns et al 1998) 
it is found that quantities such as the mean effective optical 
depth are indeed sensitive to resolution. In the SC model 
computations, although the smoothing scale does not enter 
directly as a parameter, the quantities which do enter, 7 
and erf, are dependent on it. A smaller Jean's scale will yield 
higher values of a\ , and as we can see from the N-body only 
tests of Figs. and [l|, the accuracy of our approximations 
can vary widely. The first three moments of the flux are re- 
covered within ~ 10 — 20%, for a\ ~ 4 and below (for which 
the non-linear variance in the density field, <r 2 ~ 10). The 
kurtosis of the flux has a much larger error, although as we 
will see later some non-Gaussian models have such a differ- 
ent 64, that it should still be detectable. A more direct test 
involves consideration of results from the higher resolution 
version of the simulations from set (a). This simulation has 
a higher mass resolution by a factor of ~ 12. When we use 
the same filter size as for set (a), we find that the flux mo- 
ments do not change. If we decrease the filter size by 12 s to 
0.13 /i -1 Mpc, (at which scale a\ — 3.8), the SC predictions 
for the mean flux and variance are still accurate to better 
than 10% , but the prediction for 5*3 (<f>) is —2.3, when the 



Figure 11. The mean flux (top panel), variance (middle), S3 
(triangles), and S4 (squares) for simulation set (b), as a function 
of cr 2 . We have used A = 0.8 and (3 = 1.6 to generate the spectra. 
The corresponding predictions of the spherical collapse model are 
shown as lines. All results are in real space. 




Figure 12. The mean flux (top panel), variance (middle), S3 
(triangles), and S4 (squares) for simulation set (c), as a function 
of cr 2 . We have used A = 0.8 and (3 = 1.6 to generate the spectra. 
The corresponding predictions of the spherical collapse model are 
shown as lines. All results are in real space. 



simulation value is —3.0, and the SC Si{(j)) is —1.0 when the 
simulation gives +2.4. 

One way of looking at these statistics is as constraints 
on the unknown parameters which describe the density field, 
7 and a\ . If we return to Fig. [ji] we note that the results for 
cr| are quite similar for both sets of simulations. This can be 
also seen in Fig. ^, if we look at the results for different values 
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P and A. The mean flux or the skewness seem to respond 
more sensitively to gl- We can also see this by examining 
Figures [I] and^, where for small ox the skewness of the flux 
is a much better indicator of 7 than the mean flux. This then 
changes at a\ ~ 2 where there is a degeneracy in 5*3 (4>) 
for different values of 7, while (<f>) seems to give different 
predictions. It is interesting that, although one model has 
a larger amplitude of mass fluctuations, the variance in the 
flux is not systematically higher or lower in one model than 
the other, but instead which is higher depends on A. We also 
note that the kurtosis becomes quite large for small values of 
A. Although the accuracy of our predictions for 1S4 decreases 
as erf, becomes large (see above), when of, is moderate, the 
rapid increase in S4 seen for small values of A is reproduced, 
even for values as large as 5*4 ~ 100. 

5.5 Comparisons with the mean flux held fixed 

We have seen in Fig. ^ that if we choose a specific value of 
A, the flux statistics can change fairly drastically if we add 
or remove the effects of peculiar velocities or thermal broad- 
ening. When working with observational data, the value of 
A is at best only known from estimates of the individual 
parameters in equation (§), Q b , F, H(z) to within a factor 
of 2. It would therefore be a good idea if we could fix A 
directly using Lya observations. In principle, when working 
within the formalism we have adopted in this paper, there 
are 4 unknown quantities, 7, o;, A and j3. It has already 
been found in Croft et al. (1998a) that a convenient way of 
effectively determining the correct value of A to use in nu- 
merical simulations is to choose the value which yields the 
observed mean flux (<f>). We carry out the same procedure 
here, so that when evaluating our predictions we make sure 
that {(f>) is a fixed value. In Fig. [13j we show results for three 
different values of (<f>), 0.6, 0.7, and 0.8, which are in the 
range measured from observations at z = 2 — 4. In the top 
panel, we show the value of A required for each value of {(f>) , 
as a function of of . All results are in redshift-space, except 
for ((f)) = 0.7, which we also show in real-space. We can see 
that once the mean flux is held fixed, the one-point moments 
of the flux change only little with the addition of redshift 
distortions. In particular we find that the normalized hier- 
archical moments, 53 (<j>) and S4,(<p) are practically identical 
in real and redshift-space. 

In Fig. ^ we concentrate on {(j>) = 0.7 and show the 
variation with 7, again in redshift-space. We can see that 
when the fluctuations in the underlying density field are 
large, the variance and skewness tend to asymptotic val- 
ues. As expected, more negative values of 7, will produce 
spikier, more saturated absorption features and so give a 
larger variance and skewness. To see the maximum values 
that this tendency will produce, we can consider the fact 
that 4> can only lie between and 1. There is therefore a 
limit to the level of clustering, given a specific value of the 
(4>). This will occur if a spectrum contains only pixels which 
have either c/> = 1 or c/> = (either no absorption or satu- 
rated). If the a fraction / of the spectrum has = 1, and 
the rest, a fraction (1 — /), has <f) — 0, then (d>) = f. Using 
the definition of o 2 and 53 (4>) in equations (£(]) and (p7j), 
we find that the variance in this case is 
7 1.1. 




Figure 13. Variation of the SC predictions (variance, middle 
panel and skewness, bottom panel) with erf for a fixed mean 
flux. In the top panel we show the value of A needed to give 
mean flux (<j>) = 0.6 (dotted line), 0.7 (short-dashed line), and 
0.8 (long-dashed line). All results are in redshift-space (with ther- 
mal broadening) except for the continuous line which represents 
results for (0) = 0.7 in real-space. 



and the normalized hierarchical skewness is given by 
I .\ «</>>-!) 



S 3 (4>) = {4>) 



<</>> 



1 + 



([!/«>>] - I) 2 



(38) 



(37) 



For ((f>) = 0.7, as in Fig. |l4| the maximum possible o 2 c6 = 
0.43 and S 3 (ct>) = -1.33. 



6 DISCUSSION 

The methods we have introduced in this paper are primarily 
meant to be used as tools in the study of structure formation. 
The predictive techniques for the one-point statistics should 
be most useful when combined with information on clus- 
tering as a function of scale (two-point statistics) , which we 
explore in a separate paper. As far as applying our formalism 
to observations is concerned, the one-point statistics we have 
discussed should in principle require data which resolves the 
structure in the forest, at least at the level of the Jean's 
scale, or the thermal broadening width. This would include 
Keck HIRES data (e.g., Kirkman & Tytler 1997) or other 
data with a spectral resolution better than ~ 10—15 km s — 1 . 
Use of lower resolution data effectively involves smoothing 
along the line-of-sight, and we leave the treatment of this 
anisotropic smoothing window to future work. 

At the simplest level, one could use our predictions for 
the one-point moments to compare directly to observations. 
For example, given a measurement of the mean flux (say 
{(/}) = 0.7), and a value for slope 7, one can use Fig. [14 to 
infer the predictions of gravitional instability for ^3(0) and 
5*4 ((f>), and then check them against the observed values. 
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Figure 14. Variation of the SC prediction for the flux moments, 
with a fixed mean flux, as in Fig. [l3| except this time for three 
different values of the linear slope 7. The dotted, short-dashed 
and long-dashed lines show the predictions for 7 = 0, — 1 and —2 
respectively. All results in redshift-space, and (<j>) = 0.7 for all 
curves. 



This sort of test, while being suitable for checking whether 
predictions are generally compatible with Gaussian initial 
conditions, and gravitational instability, is unlikely to be 
useful for discriminating between popular Gaussian models, 
which have similar one point flux PDFs. For example, there 
is little difference in the behaviour of the flux moments for 
two models with different values of 7 shown in Figs, [ll] and 
|l2] (see also the small differences between fi = 1 CDM and 
n = 0.4 CDM in Rauch et al. 1997). We therefore advertise 
our one-point analytic predictions as being more suitable for 
making wide searches of parameter space in order to ascer- 
tain the broad statistical trends expected in gravitational 
instability models (for example, see below for the evolution 
of the moments with redshift). Direct comparisons with ob- 
servational data will be more fruitful when carried out with 
two point statistics. We will explore these, and how the an- 
alytic methods we have developed here can be applied to 
them, in a future paper. For the moment, several obvious 
applications of our techniques for making one-point predic- 
tions suggest themselves, and we discuss these now. We also 
discuss the accuracy of the density evolution predictions, 
and compare the present work to that of others. 



6.1 Non-Gaussian initial conditions 

There is a large parameter space of non-Gaussian PDFs to 
choose from (see Fosalba, Gaztanaga & Elizalde 1998) . Here 
we will take a conservative approach and choose a model 
with mild non-Gaussianities, with hierarchical correlations, 
i.e., constant Sj, so that the cumulants kj ~ fcj -1 They 
tend to the Gaussian result as k-z 



the dimensional scaling k 



J 



k: 



J/2 



-> more quickly than 
As an illustrative ex- 



Figure 15. A comparison of Gaussian and non-Gaussian models. 
Mean flux (0) and variance (bottom) are plotted as a function 
of the linear variance <r£ for 7 = 1, A = 0.8 and (3 = 1.6. The 
short-dashed and long-dashed lines show the SC predictions for 
a Gaussian model and a non-Gaussian PT3 model respectively. 
Symbols corrspond to the Gaussian simulations, with the squares 
representing set (a) (which has 7 ~ 1) and the triangles set (d) 
(which has 7 ~ — 1.5) . The continuous lines are the perturbative 
predictions of Section 4.3 All results are in real-space. 



known chi-squared distribution (see e.g., Fosalba, Gaztanaga 
& Elizalde 1998 for details), also known as Pearson's Type 
III (PT3) PDF or Gamma PDF: 



P(p) = 



F(l/cr 2 )( C r 2 ) 1 ^ 



■ exp 



(-£)■ 



(39) 



ample, we will show results for a PDF based on the well 



for p = 1 + 8 > 0. This PDF has S3 = 2 and S 4 = 6 (in gen- 
eral Sj = [J — 1]!). The number of the chi-square degrees of 
freedom, N, in the discrete version of this distribution would 
correspond to = 2/ a 2 . It is not difficult to build a non- 
Gaussian distribution with arbitrary values of S3 or 1S4, but 
to keep the discussion simple we will just concentrate on the 
PT3 model. This model is not only mildly non-Gaussian, in 
the sense of being hierarchical, but also has moderate values 
of Sj (e.g., for comparison gravity produces S3 = 34/7 from 
Gaussian initial conditions and unsmoothed fluctuations). 

A possible motivation for introduction of the PT3 could 
be the isocurvature CDM cosmogony presented recently by 
Peebles (1998), which has as initial conditions a one-point 
unsmoothed PDF given by a chi-square distribution with 
N — 1 (i.e. a PT3 with a 2 = 2). Smoothing would in- 
troduce higher levels of non-Gaussianities through the two- 
point function, so this is a conservative approach. A more 
generic motivation follows from the arguments presented in 
Fosalba, Gaztanaga & Elizalde (1998), who find that this 
distribution plays a central role in non-Gaussian distribu- 
tions that arise from combinations of a Gaussian variables. 

In Figures |l5| - |l^ we compare the Gaussian and non- 
Gaussian PT3 predictions. As can be seen, even in this 
mildly non- Gaussian model the predictions are quite dif- 
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Figure 16. A comparison of Gaussian and non-Gaussian models. 
Wc plot the skewness 53 and kurtosis S4. The symbols and line 
types are the same as in Fig. ^ The continous lines are the per- 
turbative predictions, which are different for the Gaussian (thin 
line) and non-Gaussian (thick line) models. 



ferent and can be clearly distinguished from the Gaussian 
simulations (symbols). The SC predictions for Gaussian ini- 
tial conditions have been evaluated for 7 = — 1. The squares 
correspond to simulation set (a), which also have 7 ~ — 1. 
The triangles are from set (c) and have 7 = —1.5. This value 
of 7 is different from that used for the SC Gaussian predic- 
tions, although the simulation point for o^ is still closer to 
them than the non-Gaussian predictions. 

The leading order perturbative predictions, shown as 
continuous lines in the Figures, are the same for (</>} and o^, 
because the non-Gaussian model is hierarchical. However for 
S 3 and £4 the predictions are different. It is easy to show 
that to leading order the (hierarchical) non-Gaussian model 
predictions are given by equation (B2J) with: 



53 = S 3 (IC) + S 3 (G) 

5 4 = S4(IC) + S4(G)+4S 3 (G)S 3 (IC) 



(40) 



where S 3 (IC) and S^/C) are the initial conditions 
[Sj(IC) = (J — 1)! for PT3] and S 3 (G) and St(G) are the 
leading order gravitational values: 83(G) = 34/7 + 7 and 
S 4 (G) = 60712/1323 + 62/37 + 7/37 2 . As can be seen in 
Fig. nq the perturbative values are only reached asymptot- 
ically as <tl — > 0. The skewness S3 is significantly lower in 
the PT3 model and the kurtosis S4 is much larger, changing 
from £4 ~ —2 in the Gaussian model to S4 ~ 20 in the 
non- Gaussian one. 

Fig. [It], shows a comparion between the Gaussian and 
PT3 models, for a fixed mean flux. We noted previously (e.g., 
Fig |l^) the useful fact that there is little difference between 
the predictions with or without redshift distortions when the 
mean flux is fixed in this way. Fortunately, this is not the 
case for non-Gaussianities, as we can see here. The results 
change considerably if we assume different initial conditions. 



Figure 17. Comparison between non-Gaussian and Gaussian 
predictions when (<f>) is held fixed. We plot the optical depth 
amplitude A (top), variance ct^ and skewness S 3 (4>) (bottom) as 
a function of the linear variance a 2 L for a mean flux {<j>) = 0.7, 
7 = 1, and j3 = 1.6. The short-dashed and continuous lines show 
the predictions of the Gaussian and non-Gaussian PT3 models. 
These results are in redshift-space. 



6.2 Redshift evolution 

We can see from equation (^|) that A will change strongly 
with redshift, and that there will be a dependence on cos- 
mology through the variation in H(z). We can investigate 
how the one-point flux statistics will vary using our formal- 
ism. Wc compute the evolution of H(z) and of and show 
the results for firj = 1 and three different values of 7 in 
Fig. [l^. We can see that both o^ and S 3 ((j>) become much 
smaller as the mean absorption falls. In these figures we are 
assuming that 7 is not changing as a function of redshift, or 
equivalently that the smoothing scale is fixed in comoving 
ft -1 Mpc. 

In Fig. ^| we plot the results for three different back- 
ground cosmologies, two flat models, one with £1=1 and 
one with firj = 0.3 and Qa = 0.7, as well as an open model 
with Qo = 0.3. All models have the same 7, o; and the same 
mean flux ({(jj) = 0.66) at z=3. We can see that there is vir- 
tually no visible dependence on cosmology in the plot. The 
main source of variation in A which we have not included, 
is the evolution of the photoionization rate F, which will 
change as the population of sources for the ionizing back- 
ground changes. From Fig. [l9| it is evident that inferences 
about the evolution of the UV background should be fairly 
insensitive to the assumed cosmology. We can see why this 
is so by looking at Fig. ^(j, where we plot the evolution of A 
and of separately. If both quantities are fixed in the middle 
of the z range as we have done, then the changes over the 
range of validity of the FGPA (z > 2) are small. 
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Figure 18. The variation of the one-point statistics with red- 
shift, for three different values of 7, (dotted line), -1 (short- 
dashed line), and -2 (long-dashed line). The parameter A = 1.2 
at z = 3 in all cases, and er? = 2. Results are for £7 = 1 and are 
in redshift-space with thermal broadening. 



Figure 19. The variation of the one-point statistics with red- 
shift, for three different cosmologies. We have plotted results for 
f2 = 1 (solid line), f2o = 0.3, Ha = (short-dashed line), and 
Qo = 0.3, Ca = 0.7 (long-dashed line). In all cases, <r^ = 2 and 
(<j>) = 0.66 at z = 3. Results are in redshift-space with thermal 
broadening. 



6.3 The bias between flux and mass fluctuations 

In analogy with galaxy bias, we can define the bias of the 
flux with respect to mass fluctuations as 



(41) 



Unlike the case with galaxy bias, we can easily predict this 
quantity analytically using our formalism. We can choose 
between two sorts of bias, either the bias between the linear 
p, or the nonlinear p. In Fig. ^ we have plotted both of 
these as a function of the variance in the flux, a^. As we 
are dealing with one-point statistics in this paper, we do 
not discuss the scale dependence of bias. However we will 
do so in Paper II in this series (Gaztanaga & Croft 1999). 
We should point out though that in Croft et al. (1998a), it 
was found that the shape of a two-point clustering statistic 
of the flux (in that case the power spectrum) follows well 
that of the linear mass. The bias between the two was found 
in that paper by using a procedure which involved running 
numerical simulations set up with the power spectrum shape 
measured from observations and comparing the clustering 
level in simulated spectra with the observed clustering. In 
this paper, we can find the bias level in a simpler fashion. 
We note that for small values of the fluctuation amplitude, 
6 tends towards the values predicted by perturbation theory 
(Section 4.3). For larger values, such as those likely to be 
encountered in observations, a fully non-linear treatment, 
such as the one presented here, is needed. 




Figure 20. Variation with redshift of A and <r ; . Results have 
been plotted for 3 differemt cosmologies, Q = 1 (solid line), f2o = 
0.3, Ha = (short-dashed line), and Ho = 0.3, f^A = 0.7 (long- 
dashed line). 



6.4 Accuracy of the approximations for density 
evolution 

We have seen (e.g., Fig. [?]) that the predictions for the PDF 
have the right qualitative features (as a function of 7 and 
<tl) but do not reproduce it in all its details, even around 
5 = 0, because the SC is just a local approximation and 
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Figure 21. The bias between flux and mass fluctuations. We 
show the bias (see equation between the flux and the linear 

mass as a dotted line and the bias between the flux and the 
non-linear mass as a solid line. Both these quantities are 

plotted as a function of the observable cr^. The statistics of <f> 
have been computed in redshift space with thermal broadening, 
for {(j>) = 0.7, <x| = 2 and 7 = -1. 



does not include shear. The results of FG98 indicate that 
the statistics of the weakly non-linear density moments are 
dominated by the local dynamical contribution to the evo- 
lution of the PDF, and shear forces are subdominant (they 
tend to cancel out when taking the mean). We find here that 
a similar cancellation occurs when considering the PDF of 
the flux, (j>, even when ax ^ 1. 

Regarding the predictions for the 1-point moments 
of the flux, we have checked that the Spherical Collapse 
(SC) model yields almost identical results to the Gener- 
alized Zel'dovich Approximation (GZA), in Eq[AlC], with 
a = 21/13. This is true both in real and redshift-space, and 
also holds for the prediction of the velocity divergence 9. 
This is an interesting result because although the SC model 
is better motivated from the theoretical point of view, the 
GZA model is much simpler to implement. In particular the 
GZA model provide s us with analytical expressions for the 
PDF (i.e., equation ]A2(| and Fig. |^), which can be used in 
practice to make the predictions. 

As shown in Fosalba & Gaztanaga (1998b) the SC ap- 
proach to modelling non-linearities does not work as well 
for 6 as for p. In particular, it was found that the next to 
leading order (or loop) non-linear corrections are not as ac- 
curate, indicating that tidal effects are more important for 
6. This could partially explain why the redshift distortion 
modelling (see Fig. |^) requires the addition of an extra ve- 
locity dispersion in order to match the results of simulations. 



6.5 Comparison to other work 

The first attempts to constrain cosmology using the Lya 
forest focussed on comparisons between simulated data gen- 



erated with specific cosmological models and observational 
data, using traditional line statistics. When the simulated 
and observed spectra are decomposed into a superposition 
of Voigt-profile lines, the distribution of column densities 
of these lines and the distribution of their widths ("b- 
parameters" ) can be reasonably well reproduced by gravita- 
tional instability models (e.g., Dave et al. 1997, although see 
Bryan et al. 1998). In the context of these models, the line 
parameters do not have a direct physical meaning, as these 
statistics were intended to describe discrete thermally broad- 
ened lines. It is possible to use these traditional statistics to 
characterise the amount of small scale power in the underly- 
ing density field, for example (see e.g., Hui, Gnedin & Zhang, 
Gnedin 1998). However, statistics which are more attuned 
to the continous nature of the flux distribution and the un- 
derlying density field have advantages, as well as promising 
to be more sensitive discriminants, continous flux statistics 
can be designed to be less affected by noise and choice of 
technique than profile fitting. 

As the modern view of the Lya is essentially an out- 
growth of structure formation theory, it makes sense to bor- 
row statistical analysis techniques used in the study of the 
galaxy distribution. Unlike the galaxy distribution, however, 
the Lya forest offers a truly continuous distribution of flux, 
with no shot noise (albeit in 1 dimension) , and a well moti- 
vated theoretical relationship between the observed flux and 
the underlying mass. 

So far, analysis of spectra using such continuous flux 
statistics has mainly involved specific cosmological models, 
and direct comparison of simulations with observations. The 
mean flux, (<j>) is the most obvious flux statistic to calculate. 
Its measurement from observations has been carried out by 
several authors (e.g, Press, Rybicki & Scheider, 1993, Zuo & 
Lu 1993), and there is an extensive discussion in the litera- 
ture about what is usually quoted as the mean flux decre- 
ment, Da = 1 — (0), or the mean effective optical depth, 
"off = —ln(4>) 

The probability distribution of the flux has been inves- 
tigated by Miralda-Escude et al. (1996), Croft et al. (1996), 
Cen (1997), Rauch et al. (1997), Zhang et al. (1998), and 
Weinberg et al. (1998). Other statistics such as the two point 
correlation function of the flux have been introduced (Zuo 
& Bond 1994), the power spectrum of the flux (Croft et al. 
1998a, Hui 1999, the two point pdf of the flux (Miralda- 
Escude et al. 1997), and the number of times a spectrum 
crosses a threshold per unit length (Miralda-Escude et al. 
1996, Croft et al. 1996, Weinberg et al. 1998). Methods have 
been developed to reconstruct properties of the underlying 
mass distribution, such as the matter power spectrum (Croft 
et al. 1998a, Hui 1999, using our theoretical assumptions for 
the relation between mass and flux. A technique for carrying 
out a direct inversion from the flux to the mass distribu- 
tion has been described by Nusser and Haehnelt (1998). In 
the present paper, we emphasise the use of statistics which 
have been used extensively in the study of galaxy clustering, 
in particular the higher order moments (e.g., Gaztanaga & 
Frieman 1994). These statistics, and their behaviour when 
used to quantify the evolution of density perturbations in 
the quasil-linear regime have been the subject of much at- 
tention. It would seem that extending their use to the study 
of the Lya forest may offer us a good chance to combine our 
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knowledge of gravitional instability with that of the IGM 
and in doing so make progress in both disciplines. 

On the predictions side, many pieces of analytic work 
have been carried out which incorporate the dominant phys- 
ical processes involved in producing high-redshift Lya ab- 
sorption (processes summarized in Section 2). The studies' 
most important differences have been in the schemes used 
to follow the evolution of density perturbations. These have 
included linear theory (Bi 1993, Bi, Ge & Fang 1995) the 
lognormal approximation (Gnedin & Hui 1996, Bi & David- 
sen 1997), and the Zel'dovich Approximation (McGill 1990, 
Reisennegger & Miralda-Escude 1995, Hui, Gnedin & Zhang 
1997). Unlike these approximations, the SC model used in 
this paper is able to reproduce exactly the perturbation the- 
ory results for clustering. This accuracy makes it useful for 
calculating high-order statistics of the flux, in our search for 
the signatures of gravitational instability. It is important to 
realize that we have not used the SC model to make simu- 
lated QSO Lya spectra, but that we have used its predic- 
tions for the properties of the mass to predict the statistics of 
the Lya flux. With such an approach (similar to that taken 
by Reisennegger & Miralda-Escude 1995) we can quickly 
and easily vary parameters in order to explore for example 
the dependence of a particular statistic on redshift (Section 
6.2). Some tasks which previously required numerical sim- 
ulations, such as finding the bias between density and flux 
fluctuations (e.g., Croft et al. 1998ab), can be carried out 
analytically. 

The fully non-linear analysis we have described in this 
paper will allow one to carry out many analyses of cluster- 
ing where the precise relationship between the statistics of 
the mass and the flux is important. This includes attempts 
to constain the cosmic geometry from the clustering mea- 
sured between adjacent QSO lines of sight (e.g., McDonald 
& Miralda-Escude 1999 Hui, Stebbins & Buries 1999. In such 
situations, a non-linear theory of redshift distortions in the 
Lya forest and of the bias between the fluctuations in the 
observed flux and the mass, such as we have presented in 
this paper should be very useful. 



7 SUMMARY AND CONCLUSIONS 

We have presented a fully non-linear analytical treatment of 
the one-point clustering properties of the high-z Lya forest 
in the gravitational instability scenario. The formalism we 
have presented should prove to be a useful tool for studying 
the forest, and has immediate application to the calculation 
of two-point statistics (see Paper II). The two main ingredi- 
ents we have used are the Spherical Collapse model (SC) or 
shear-free approximation for the evolution of density pertur- 
bations, and the Fluctuating Gunn-Peterson Approximation 
for the relation between density and Lya optical depth. The 
predictions for the one-point clustering of the mass made us- 
ing the SC model depend only on two parameters, a L and 7. 
These are, respectively, the linear variance of density fluctu- 
ations, and the local slope of the linear correlation function. 
In the FGPA, the relation between the mass distribution 
and Lya forest optical depth is largely governed by one pa- 
rameter, A, which can be set by appealing to observational 
measurements of the mean flux. The predictions of the SC 
model for the density are typically quite non-linear (a 2 ~ 2 



or more). While these predictions are not expected to be 
accurate for the high density tail of the distribution, the 
weighting of the FGPA relation means that the statistics of 
the flux are governed by the (quasi-linear) density regime 
where the SC is accurate. The Lya forest is therefore well 
suited to study using such an approximate analytical tech- 
nique. 

We note that the analytical predictions can be used in 
tests of the picture of Lya forest formation and the appli- 
cability of the FGPA. With the extra information afforded 
by the two-point statistics and considering the evolution of 
clustering as a function of redshift, it will be possible to 
look for the signatures of any deviation from the theroreti- 
cal picture. Consistency tests for the gravitational instability 
scenario include checking the evolution of the moments as a 
function of redshift, the scaling of the hierarchical moments, 
and their dependence on scale. 

We plan to use our predictive techniques in future work 
to extract information from observations, using both one- 
point and two-point statistics. In the present paper, we 
have concentrated on developing an analytical framework for 
studying the clustering of the transmitted flux in the forest 
region of QSO spectra. Some of the results of our present 
exploration of one-point statistics include the following: 

• Using our formalism we are able to estimate the bias be- 
tween mass and flux fluctuations without resorting to simu- 
lations. 

• We can make predictions for the clustering properties of 
the Lya forest flux in both Gaussian and non-Gaussian mod- 
els. We find large differences between the two in an example 
case. 

• In the limit of small fluctuations, our non-linear analyt- 
ical treatment converges to the same results as those from 
Perturbation Theory calculations. 

• For larger fluctuations, where Perturbation Theory is no 
longer valid, we find our treatment to give accurate re- 
sults compared to statistics evaluated from N-body simu- 
lated spectra. These predictions are most accurate when the 
linear variance of the density field is ~ 4 and below. For val- 
ues above this, the qualitative behaviour of the high order 
moments is reproduced. 

• We can follow the evolution of the one-point statistics of 
the flux as a function of redshift. We find that the differ- 
ence between predictions for different cosmologies is small, 
so that comparison with observations should be useful in 
constraining the evolution of the ionizing background inten- 
sity. 

• If we normalise our predictions so that the mean flux is 
held fixed (for example to the observed value) , we find that 
the statistics of the flux are relatively insensitive to the ef- 
fects of redshift distortions induced by peculiar velocities or 
thermal broadening. This is most valid for the higher order 
normalised hierarchical moments. 

The high-z Lya forest is amenable to study using an- 
alytical treatments, a fact which gives it great value as a 
probe of structure formation. Application of analytical the- 
ory to measurements from the many observational datasets 
currently available promises to reveal much, both about the 
validity of our theoretical assumptions, and about cosmol- 
ogy- 
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APPENDIX A: NON-LINEAR MAPPING 
RELATIONS 

Al Unsmoothed relations 

A 1.1 The density 

For small linear fluctuations, a generic mapping can be ex- 
pressed in a Taylor series: 



(Al) 



where Uk are the coefficients that define the transformation 
in this limit. The above mapping can be used to predict the 
PDF of the evolved field, or directly predict its cumulants. 
For example, for the skewness we have: 

(5 3 } c = {(G - {G}) 3 } = 3u 2 ai + O(at) (A2) 

where we have made use of the Gaussianity of Sl ■ In general 
we have: 



<<5 J ) C = Sj a 2 L (J - 1] +0(a 2 L J ). 



(A3) 



For the skewness S3 — 3^2, and in general Sj can be 
given in terms of Uk (e.g., see FG98 or the original work 
by Bernardeau 1992, which is in terms of the generating 
functionals) . 

For the SC (see below) we have: 

34 
21 " 
682 

189 

446440 



V2 = 



1.62 
- 3.61 



V4, = 



43659 
8546480 
243243 



10.22 
- 35.13 



(A4) 



These numbers give the correct leading PT contribution 
to the Jth-order cumulants {5 J )c for Gaussian IC, e.g., 



S3 = 34/7 (Peebles 1980, §42). The important point to no- 
tice here is that although the local mapping is not the exact 
solution to the evolution of S (which is in general non-local) , 
it gives the correct clustering properties in the weakly non- 
linear regime. This is because of the symmetry involved in 
taking the ensemble average (...) (Bernardeau 1992). This is 
also a good approximation for next to leading terms (FG98) 
and also for non-Gaussian initial conditions (Gaztanaga & 
Fosalba 1998). 



A 1.2 Velocity divergence 

The density mapping can be used to predict the velocity 
divergence, defined here as: 



= -V-v 



(A5) 



where v is the peculiar velocity field. We can use the conti- 
nuity equation: 

dp 



dt 



Hp6 = 



(A6) 



to express 9 as a Lagrangian mapping of Sl' 

,A7) 

where fa = d In D/d In a, comes from applying the chain rule 
to the derivative of the linear growth factor: Sl = D(£)5ic- 



A 1.3 The SC model 

For large p the spherical collapse mapping 8 — G[8l] can 
only be expressed in a parametric form through an auxiliary 
variable ip: 

+ 9 (i/> - sin^) 2 



y 2(l-cosV) 3 ' 
Si = l^-sm^f 3 , 

for Sl > 0, which we call 5t , and 



(A8) 



9 (sinhV - VO 2 



2 (coshi/> 
_3 3 
5 l 4 



1)3' 

(sinhV - i>)] 2/3 , 



(A9) 



for Sl < 0, which we call 5^. This parametric form 
complicates the analytical predictions, but it can be imple- 
mented numerically. For ip — 2-k the transformation becomes 
singular and p and <5 diverge. This is the first collapse which 
occurs at S L = 3/5(3tt/2) 2/3 ~ 1.686, as illustrated in Fig. 
Al. It then bounces back and recollapses in a periodic fash- 



A1.4 The GZA model 

Another model we consider for the local transformation is 
the Generalized Zel'dovich Approximation (GZA): 



1 + .5: 



1 



Sl 
a 



(A10) 
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Figure Al. Local mappings for positive (top right-hand corner) 
and negative fluctuations (bottom left-hand corner). The short- 
dashed line shows the (unsmoothed 7 = 0) spherical collapse (SC) 
result, which has a first bounce at 5; ~ 1.69. The continuous 
line (which almost covers the short-dashed line) is the GZA with 
a = 21/13 dimensions, which becomes singular at Si ~ 1.62 and 
has its first two derivatives equal to those of SC. The three long 
dashed lines below show the same GZA smoothed with increasing 
smoothing indices of 7 = —1.5, —3, —5. These curves are difficult 
to distinguish from the corresponding SC mappings smoothed 
with the same 7, which are displayed as dotted-lines. 



which is a symmetric version of the ZA in a dimensions and 
was also considered by Protogeros and Scherrer (1997). The 
case a = 3/2 was introduced by Bernardeau (1994) as a good 
approximation to the SC model in the weakly non- linear 
regime. This value of a gives vi ~ 1.67 in equation (Al), 
compared to the exact value V2 = 34/21 ~ 1.62. Here we will 
concentrate instead on the case a — 21/13, as it reproduces 
exactly the SC result to second order: V2 = 34/21, and there- 
fore gives the exact skewness, 5*3 — 34/7. This value of a also 
is a better approximation than a = 3/2 around the singular- 
ity in 5, which occurs for Sl = a — 1.62, closer to Sl — 1.689 
than Sl — a — 1.5. The higher order derivatives for the case 



a = 21/13 are vz — 3.62, ^4 ~ 10.35 and v$ = 35.99, which 



should be compared to the SC values in equation (A4). An- 
other model of interest is a = 3/5(3tt/2) 2/3 ~ 1.686 which 
reproduces the behaviour exactly around the singularity, and 
gives V2 — 1.59, which is closer to the PT value in the SC 
than a = 3/2. 



The velocity divergence, equation (A7), in the GZA has 
the simple form: 



= -fa Sl P 1/a = -fa 6 



l- S -± 

a 



(All) 



An important problem with these relations is that there 
is more than one value of Sl for a given value of S. This might 
be important when it comes to obtaining the PDF of 8. In 
practice this branching problem can be solved by assigning 
to 5 the probability associated with the different branches in 
Sl- However, we have checked that this is not important for 



smoothed fields and this allows us to ignore this branching 
problem from now on. 



A2 Smoothing effects 

The above local relations correspond to unsmoothed fluc- 
tuations. We will focus instead on fields smoothed with a 
top-hat filter, defined as, 



Wth(X,R) 



1 if \X\<R, 



and zero otherwise, where R is the smoothing radius. The 
volume corresponding to radius R is V — 4-k/3R 3 . For top 
hat smoothing, the statistical properties of the smoothed 
fields can be obtained using the following prescription (see 
also Bernardeau 1994, FG98). Consider an initial (La- 
grangian) fluctuation of infinitesimal mass mo extending 
over some volume Vo = 47r/3i?o- The statistics of the IC 
must be given as an input to estimate the corresponding 
evolved, non-linear, values after gravitational growth. This 
is also true for the smoothed case, so that the shape and 
amplitude of the initial variance must be given or set 'by 
hand'. Consider first the case were we want the initial fluc- 
tuations smoothed within a radius R, to have a power-law 
variance: 



~ 2 2 
a = <7 



(ro) 



(A12) 



where Ro and erg relates to the "unsmoothed" amplitude. In 
the linear regime, this sets the amplitude of linear fluctua- 
tions: 



R V' 2 



(A13) 



where D(i) is the gravitational growth factor and <5o is some 
(unsmoothed) initial seed whose variance is a 2 , so that in 
the limit R — > Ro we have a\ — > D 2 a 2 = a\. 

On the other hand, in the local picture of evolution, a 
given fluid element mo is isolated, so that for a mean density 
n, the smoothed overdensity is: 



m / Ro Y 

p= ^v = \.-r) ' 



(A14) 



as the IC are perfectly homogeneous (mo = nVo), and the 
mean density n does not change with time. 

Puting the above equations together we have: 



P = 



6/7 



(A15) 



which shows that smoothing acts like a (implicit) La- 
grangian mapping. Thus, given the unsmoothed mapping 
p ~ G(Sl), we can find the corresponding smoothed map- 
ping by solving the implicit relation: 



g[p 1/6 s l 



(Ai6) 



where 7 is the slope of the initial or linear smoothed vari- 
ance: 

7=^4- (A17) 
dlog R 

This can be easily generalized to non power-law relations 
and reproduces the original result of Bernardeau (1994). 
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The velocity divergence can be thought as a different 
Sl mapping, i.e., equation (A7). T hus t he smoothed results 



can be obtained as with equation (A16) 



7/6-1 



6 = Q v [8lP 



A2.1 Smoothing in the SC model 



(A18) 



In principle, the smoothing in p is difficult to implement an- 
alytically for the SC , because Q is given by the implicit rela- 
tions Eqs. A8|-[A9|. However, the smoothing can be carried 
out easily in terms of the smoothed Sl, sim ply by noticing 
that from the above relation (equation ]A1E[ ) we have: 

(A19) 



where ip is such that p = p[ip]. The smoothed PDF of p will 
be given in terms of the smoothed PDF of the IC, which for a 
Gaussian field is also a Gaussian. So we just have, following 
Eq§: 



P(P) 



Pic{5l 



dS L 
dp 



Pic(5 L [j>] 

pm 



dS L 

dip 



dp 
Thp 



, (A20) 



with 8h[ip\ given by equation ( A19 ) . Ca re must be taken 
to use the correct relation, equation (A8) or equation ( A9), 
depending on the sign of Sl- For example, to estimate the 
moments of p we have: 



(P J ) 



P(p)p J dp 



+ 



1 

77 
1 

77 



Pic (Sj 



P + \% 
Pic 01 



dsi 



dip 
d57 



dip + (A21) 



dip 



[4>} J dip, 



where is the normalization factor, so that (1) = 1. The 
superscript + denotes our use of equation (A8) and _ our 
use of equation ( A£ ) . 



For example, for the evolved (smoothed) one-point PDF as 
a function of p, we have: 



P(P) 



1 

77 



exp 



l/a 



/(2a 2 



(A25) 



-i/a-7/6-2 



7 r l/a 
6 a[p 



for Gaussian inital conditions. These sort of compact rela- 
tions were previously described by Protogeros & Scherrer 
(1997). Fig. |i] shows a comparison of the above PDF to the 
PDF from simulations and from the SC model. 

The velocity sm oothi ng can be obtained from equation 



(A18) and equation (All) 



= -fa 5lP 



i/a+7/6 



= -fa Sl 



1-^ 



(A26) 



Note that although the SC and the GZA give a good 
approximation to gravitational dynamics for small S, it is 
likely and expected that these approximations break at some 
stage for large fluctuations. As mentioned before, this break 
is not likely to be important for our present purposes, as 
the statistics of the absorption features in the Lya forest 
are dominated by the small fluctuations. 



A2.2 Smoothing in the GZA model 



From now on, unless stated otherwise, we use S and p to refer 
to the smoothed fields. The unsmoothed case corresponds to 
7 = 0. Th e sm oothed GZA mapping is given implicit ly by 
equation (A16), which for the GZA case of equation (A10) 
yields: 



P - 



Sl p 



.1/6 



(A22) 



which nicely maps Sl G [—00,00] into p £ [0, 00]. The re- 
sulting smoothed mappi ng f or the GZA for a — 21/13 and 
7 = —2 is shown in Fig. Al. 



It turns out that for the GZA, even after smoothing, we 
can give analytical expressions for the PDF simply by using 
Eq||l with: 



1 



dS L 
dS 



o7 /6+l/c 



-l/a-7/6 _ 7 £ 

a U L 



(A23) 
(A24) 
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